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■ Abstract 

We consider the (2 + 1) flavor Polyakov Quark Meson Model and study the effect of including fermion vacuum 
fluctuations on the thermodynamics and phase diagram. The resulting model predictions are compared to the 
recent QCD lattice simulations by the HotQCD and Wuppertal-Budapest collaborations. The variation of the 
thermodynamic quantities across the phase transition region becomes smoother. This results in better agreement 
with the lattice data. Depending on the value of the mass of the sigma meson, including the vacuum term results 
in either pushing the critical end point into higher values of the chemical potential or excluding the possibility of 
a critical end point altogether. 
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1 Introduction 

I With the advent of ultrarclativistic heavy ion experiments such as RHIC (BNL), LHC (CERN) and the planned 
■ future CBM experiment at the FAIR facility in Darmstadt, research on the properties of strongly interacting matter 
has entered into a fascinating era. Quantum Chromodynamics (QCD) which is the theory of strong interactions 
predicts that at high temperature and baryon density normal hadronic matter undergoes phase transition to a chirally 
symmetric phase of Quark Gluon Plasma (QGP) [1-5]. 
t-H ■ QCD at this temperature and density is a strongly coupled theory and hence one cannot use perturbative methods 
^ | to study the details of this phase transition. Non-perturbative methods like first principle Lattice QCD (LQCD) Monte 
• ■ Carlo simulations give us important insights into various aspects of the phase transition. However LQCD suffers from 
| the notorious sign problem at non zero baryon density [6,7]. Although several methods have been developed [6,7] 
to bypass the sign problem at small baryon chemical potential /ib, a satisfactory solution to the sign problem for all 
values of [Ib still eludes us. An alternate approach is to study various phenomenological models whose phase diagrams 
possess the essential features of QCD. These models serve to complement LQCD calculations and may provide us with 
an intuitive understanding about regions of the phase diagram currently inaccessible to LQCD. 

The QCD Lagrangian for Nf flavors of massless quarks has the global symmetry SC/l (iV/) x SUn(Nf) which 
spontaneously breaks into SUy (Nf) in the low energy hadronic vacuum by the formation of chiral condensate together 
with Nf — 1 massless Goldstone bosons. In the opposite limit of infinite quark mass, QCD becomes a pure SU (3) 
gauge theory. The low energy vacuum that possesses a center symmetry Z (3) under the color gauge group (confined 
phase) gets spontaneously broken in the high temperature/baryon density regime (deconfined phase). The Polyakov 
loop which is vanishing in the confining phase and becomes non zero in the symmetry broken phase serves as the order 
parameter of the confinement-deconfinement phase transition [8] . In real life we have dynamical quarks with non zero 
masses and hence neither of the symmetries discussed above arc exact. Thus the chiral condensate and the Polyakov 
loop only serve as approximate order parameters for the chiral phase transition and confinement-deconfinement phase 
transition respectively. 
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There have been numerous attempts to study the QCD phase diagram using effective chiral models like the 
Nambu-Jona-Lasinio (NJL) [9-12], linear sigma [13-15] and the quark-meson (QM) [16] models. These models possess 
QCD like chiral symmetry breaking pattern. Later, variants of these were developed to incorporate the confinement- 
deconfinement transition, namely PNJL [17-20], PLSM and PQM [21,22] models. There is a basic difference in the 
mechanism of chiral symmetry breaking in the QM/PQM models compared to NJL/PNJL models. In the latter 
case, the fcrmion vacuum term which is the contribution due to the infinite Dirac fcrmion sea causes the symmetry 
breaking. In QM/PQM models, the spontaneous symmetry breaking is generated by the mesonic potential and the 
fcrmion vacuum term is usually neglected. Recently in the 2 flavor PQM model it has been shown that the order of 
the phase transition in the massless limit changes from first order to crossover on adding the vacuum term [23] which 
also modifies the shape of the isentropic trajectories near the critical end point (CEP) [24]. In yet another study in 
the presence of a magnetic field [25], the phase diagram is shown to be considerably affected by the vacuum term. 
More recently, the effect of the vacuum term in the 2 flavor case on the thermodynamics as well as the phase diagram 
was extensively studied [26]. 

In this paper we will consider the PQM model together with the vacuum term (PQMVT) in the (2 + 1) flavor case. 
We first compute the thermodynamic potential in the presence of the vacuum term. Then we calculate thermodynamic 
quantitites like energy density e, pressure p, entropy s, conformal symmetry breaking measure A, speed of sound c 2 s 
and specific heat capacity cy . We also compute several second order quark number susceptibilities (QNS). In each case 
we compare our results with those obtained from LQCD and comment on the modification in the model prediction 
due to the addition of the vacuum term. We also present the phase diagram of PQMVT in the T — (1b plane. The 
organisation of the paper is as follows: In Sec. 2, we discuss the details of the PQMVT model and its parameters. In 
Sec. 3 we present our results and elucidate the effect of the vacuum term on thermodynamics and phase diagram of 
PQMVT. Finally, we summarise and conclude in Sec. 4. 



2 The Model 

We will be working in the generalized three flavor quark meson linear sigma model which has been combined with 
the Polyakov loop potential [22,27,28]. In this model, quarks denoted by qf — (u,d,s) T come in Nf = 3 flavor and 
N c = 3 color degrees of freedom and are coupled to the SUl(3) x SUr(3) symmetric mesonic fields M as well as a 
spatially constant temporal gauge field A 11 . Our starting point is the following model Lagrangian 

C = Tr (d^d^M) + q f {i-fDp - gM 5 )q f - U M (M) - U P ($, 3>, T) , (1) 

where nine mesons in the scalar (a a , J p = + ) and pseudoscalar (n a , J p = 0~) sectors are together put into the 3x3 
complex matrices M and M 5 given by 

M = T a ^ a = T a (a a + in a ) (2) 
M 5 - T„6a = T a (a a + i~f 5 ir a ) (3) 

where T a = 4f represent the 9 generators of U(3) with A = ^/| 1. The generators follow t/(3) algebra [T a ,Ty = 

ifabcTc and {T a ,Tb} — d a ^ c T c where f a b c and d a bc are standard antisymmetric and symmetric structure constants 

respectively with f a b = and d a b = ^/| 1 (5 ab and arc normalized as Tr(T a T b ) = 

In (1), Um (M) and Up (<&, $,T) are the mesonic and Polyakov loop potentials respectively. Polyakov loop field 
3>(x) is defined as the thermal expectation value of color trace of the Wilson loop in the temporal direction 

$(x) = ^ r (Tr c L(x)), (4) 



where L (a?) is a matrix of the SU C (3) color gauge group. 



L(x) — "Pexp 



f 

i I drA (x,T) 



(5) 



Here V denotes path ordering, Aq is the temporal vector field and (3 = T 1 [8]. Also, $ = -^-(Tr c L^ (x)) is the 
Hermitean conjugate of $. The coupling of quarks with the gauge field is implemented through the covariant derivative 
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Dp = 8^ — iAp and = S^qAq (Polyakov gauge), where A^ = g s A a ^\ a /2. \ a are the Gell-Mann matrices in the color 
space where the index a runs from 1 • • • 8 and g s is the SU C (3) gauge coupling, g is the flavor blind Yukawa coupling 
that couples the three quark flavors with the mesons a a and ir a of the meson matrix M5. 
The mesonic potential Um (M) has the following form 



U M (M ) = m 2 Tr(M t M) + Ai [Tr(M^M)] 2 

+A 2 Tr (M^M) 2 - c[det(M) + det(M^)] 
-Tr [H(M + M f )] . 



(6) 



While the 't Hooft determinant term breaks explicitly the U (1)^ symmetry, the SUl{3) x SUp(3) chiral symmetry 
is explicitly broken by H = T a h a and the £ field picks up a non zero vacuum expectation value £. In order to give 
the correct vacuum quantum numbers after chiral symmtery breakdown, ho, /13 and h$ are the only allowed non zero 
parameters in H. In this work we retain isospin symmetry which further sets /13 = 0. This leads to the 2 + 1 flavor 
symmetry breaking scenario with nonzero condensates do and &%. 

The effective Polyakov loop potential Up (<&, $,T) is constructed to reproduce the thermodynamics of pure gauge 
theory on the lattice. The simplest Z(3) symmetric polynomial form based on a Ginzburg-Landau ansatz is given 
by [19] 



W Poly (*,$, T) ^ b 2 (T) 



(7) 



T 4 2 

Later an improved ansatz was proposed [29, 30] by including the VanderMonde term which is the Jacobian of trans- 
formation from the matrix valued field L to the complex valued field $. Here we closely follow [30] 



— ^ = —j^ - k log p - 6$$ + 4 ($ J + $ J ) - 3 ($$) 



(8) 



The logarithm term constrains the $ and $ values to smaller than 1. Details of the above mentioned model parameters 
will be discussed in Section 2.2. 



2.1 Grand Potential 

The grand canonical partition function for a spatially uniform system in thermal equilibrium at finite temperature T 
and quark chemical potential /x/(/ = u, d, s) is given by 



Trexp[-/?(7} 



f—u.d,s 



= JjJ D<r a Dir a J VqVq exp - J dr J d 3 x 



f—u.d,s 



(9) 



where V is the three dimensional volume of the system, and j3 = y. In general, the quark chemical potentials 
corresponding to the three conserved flavors up(u), down(d) and strange(s) are different. In the 2 + 1 flavor setup, the 
SUv{2) symmetry is unbroken and thus the quark chemical potentials for u and d quarks become equal yu„ = [i d = \i x 
while the strange quark chemical potential is fi s = fj, y . 

We evaluate Z in the mean- field approximation [31,32]. In this approximation, thermal and quantum fluctuations 
of the meson fields are neglected and are replaced by their expectation values (M) = To(cro) + T%{o%) while quarks 
and antiquarks are retained as quantum fields. Now following the standard procedure [18,19,21,33] one can obtain 
the expression of the grand potential as a sum of a pure gauge field contribution Up ($, $, T), a meson contribution 
Um {(M)) and a quark- ant iquark contribution ttq q 



TlnZ 



= U M ((M)) +U p ($,$,T) +SI 



<?<? 



(10) 



In order to study the 2 + 1 flavor case, a more convenient basis is the nonstrange-strange (x, y) basis obtained 
from the original singlet-octet (0,8) basis by the following basis transformation, 



J_ ( V2 1 



vo 



(11) 
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applied on the scalar, pseudoscalar mesons and external fields, that is v e {a, ir, h}. In this new basis the nonstrange 
and strange quarks and antiquarks decouple and the quark masses become 



9 Y , 



TO,, 



>T2 



The mesonic potential in this nonstrange-strange basis reads, 

U(a x ,a y ) = W M ((M» 



+ |(2A 1 +2A 2 )^ , 
while the quark- ant iquark contribution is given by, 



— \v x + v y ) - n x a x - hyOy - -^-j=<j x a y 
1 



n^n = n v - + n th 

uuqq qq 1 qq 



where 



il qq 



-2N C £ / 



f—u,d,s * 



th 



-2T £ 

f—u,d,i 



d 3 P 

(2^)3 

(2^ 



In 5^ + In ^rj 



5^ and <fy are defined after taking trace over color space 

g+ = [l + 3$ e - B / + / T + 3$ e - 2£ / /T + e - 3is / /T 
g- = [l + ^e- E f' T + ^e- 2E J' T + e^ E J' T 

Here we use the notation E* = Ef =p [if where Ef is the single particle energy of a quark/antiquark. 



Ef = yp 2 + to/ 2 



(12) 



(13) 
(14) 

(15) 
(16) 



(17) 



(18) 



At zero temperature and chemical potentials, flg q gets a contribution only from fi^. This is the fermion vacuum 
contribution which is usually neglected. Wc shall study the effect of this term on the thermodynamics. Employing 
dimensional regularisation to regularise the diverging integral in (15), we obtain [23,26] 



TO/ 



f—u,d,s 



A 



(19) 



where A is the regularisation scale parameter. In Appendix A, we show that the thermodynamic potential is indepen- 
dent of A. Thus all physical observables are independent of the choice of A [23]. 

The quark condensates (<r x ), (a y ) and the Polyakov loop expectation values ($), ($) are obtained by searching for 
a global minima of the grand potential at a given value of temperature T and chemical potentials [i x and [i y , 



dn _ an _ on _ on 



= . 



<7 a! = <^>,<7 w = <<7 w >,*=(*>,4=(*> 



(20) 
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2.2 Model Parameters 



The model parameters are fixed by inputs from lattice and experiments. The parameters of the mesonic potential 
namely to 2 , Ai, A2, c, g, h x and h y are fixed by using the following experimentally known quantities as in [15,32]: 
pion mass to^ and pion decay constant /„., kaon mass ttlk and kaon decay constant fx, average squared mass of i] 
and rf mesons, (m 2 + to 2 ,) and sigma mass m a . g is fixed by using a light quark constituent mass m x — 300 MeV 
in (12). This gives a strange quark constituent mass m s w 433 MeV. It turns out that A2 is the only parameter 
that depends on A. However, this dependence cancels with that of f2™ s (A) to yield a A independent thermodynamic 
potential. Thus, A is an arbitrary parameter here and we do not need any further experimental input to fix it. In 
Appendix A we discuss in detail the fixing of these parameters. Experimentally the mass of the a meson m a lies in 
the range 400 — 1200 MeV [34] while recent studies point to m a <~ 400 — 500 MeV [35]. All our results are obtained 
with to ct = 400 MeV. In Table 1, we have given the parameters used in this work. 

We now discuss the parameters of the gauge potentials U*p \ y ($, l>, T) and Wp i y -VM($, ®, T) used in this work. 



^Poly-VM (*,*,T) 



rp4 



W Poly ($,$,T) 

7^4 --Klog 



1 - 6$$ + 4 ($ 3 + $ 3 ) - 3 ($$) 



where Up \ y is given by 
Here b 2 (T) is given by 



rp4 



a 3 



(21) 
(22) 
(23) 



The parameters of b 2 (T), 63 and 64 are obtained by fitting to the pure gauge lattice data in [19]. 



ao = 6.75 , 
a 2 = 2.625 , 



ai = 

«3 



-1.95 , 
-7.44 



63 = 0.75 
64 = 7.5 



The remaining parameters To and k are fixed by comparing the model predictions to lattice data of the full (2 + 1) 
QCD. Tq = 270 MeV is the phase transition temperature in the pure gauge theory. In the presence of dynamical 
quarks, fermionic contributions modify the running coupling of QCD [36,37] which can be accounted for by a smaller 
value of T [21,36,37]. While both HotQCD as well as WB groups agree that the chiral symmetry restoration takes 
place through a smooth crossover [38-41], there is still a disagreement on the relative values of T x and Tj, the 
transition temperatures associated with chiral and confincmcnt-dcconfinement crossovers respectively. The HotQCD 
lattice data [42,43] with which we compare our model predictions show a coincidence of both transitions. On the 
other hand, WB sees T x < Td [40,44-46]. However, recently HotQCD have reported lower value of T x [47] which is in 
agreement with that of WB. We observe that by varying To it is possible to change the relative values of T x and Td 
in the model. We work with two values of T : firstly, T = 270 MeV. In this case T x < T d . We call this ModelWB. 
Secondly, T = 210 MeV. In this case T x ~ Td- We call this ModelHotQCD. The remaining parameter k is adjusted 
to describe the lattice data of the scaled conformal symmetry breaking measure A/T 4 = (e — 3p) /T 4 . The two lattice 
groups differ in the height of the peak of A [44]: the peak height is about 50 % larger in HotQCD. Here we use 
different values of n in ModelWB and ModelHotQCD to describe the WB and HotQCD data respectively. One could 
have obtained To and k by precision fitting of the LQCD data. Here we do not undertake such a procedure as the 
exact quantitative status is going to change as more refined LQCD data become available. Here we merely focus on 
comparing the trends between model and LQCD rather than exact matching. We summarise below the values of all 
the model parameters in Table 1. 



Model 


to (7 [MeV] 


to 2 [MeV"] 


Ai 


A 2 (A) 


c 


h x [MeV 3 ] 


h y [MeV 3 ] 


9 


T [MeV] 


K 


ModelHotQCD 


400 


80647.587 


-8.165 


138.45 


4801.95 


1.785 x 10 e 


3.805 x 10 y 


6.5 


210 


0.1 


ModelWB 


400 


80647.587 


-8.165 


138.45 


4801.95 


1.785 x 10 b 


3.805 x 10 Y 


6.5 


270 


0.2 



Table 1: The parameter sets obtained with A = 200 MeV. 
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Figure 1: Plots of normalised {a x )(red), (a y ) (green) and ($)(blue) vs T/T x at \i x — \i y — as obtained in PQMVT 
with Wp iy_vM as the Polyakov potential. The dashed curves refer to ModelHotQCD while the solid curves refer to 
ModelWB. The solid black curve in each case is obtained in PQM with the ModelWB parameter set. 

3 Results 

Having obtained the grand canonical potential (T, jjt x , n y ) in the presence of the vacuum term, we shall now compute 
various thermodynamic quantities and study the effect of inclusion of the vacuum term on these quantities. For a 
given (T, n x , n y ), at the mean field level, f2 (T, /i x , fj, y ) is only a function of the condensates a x , a y , $ and $ which 
are determined by solving the gap equations (20) simultaneously. In Fig. (1) we have plotted the condensates (cr x ), 
(<r y ) and (<£>) obtained in ModelHotQCD and ModelWB at zero chemical potentials (n x = /j, y = 0) with Wp D i y _vM as 
the Polyakov potential. We do not plot ($) since at zero chemical potential ($) = ($). While the quark condensates 
come out similar in both the models, using a higher value of To and k in ModelWB results in a delayed confinement- 
deconfinement crossover as well as a suppressed value for the Polyakov condensate as compared to ModelHotQCD. 
For large T, ($) in ModelWB falls short by 8% compared to that of ModelHotQCD. In order to compare the results 
between PQM and PQMVT, we have also plotted the results as obtained in PQM with the ModelWB parameter set. 
In PQM, (a x ) drops sharply across the transition region accompanied by an unsmooth structure in (a y ) and ($). With 
the inclusion of the vacuum term, the drop in (a x ) across the crossover is much more gentle and the jagged structures 
in the transition region in (a y ) and ($) are washed away. One can compare these condensates as obtained in PQMVT 
with those computed in LQCD. 

In the left panel of Fig. 2, we have plotted the Polyakov condensate ($) obtained for both the parameter sets with 
^Poiy-VM as the Polyakov potential and compared with LQCD [43,45]. We have also plotted ($) with Up \ y as the 
Polyakov potential and ModelWB as the parameter set for comparison. Although the VanderMonde term in Wp i y _vM 
suppresses the ($) value, the model prediction is still much higher as compared to lattice. The use of a potential that 
yields first order phase transition in case of pure glue theory results in a more rapid rise of the Polyakov condensate 
not seen in lattice as also reported in PQM studies [28]. 

A direct comparison of (a x ) and (a y ) with lattice data is not possible. Unknown normalisaton factors need to be 
first removed. One such observable with the desired chiral limit and which also acts as an order parameter for chiral 
symmetry breaking is [42,43] 

A, s (T) = X J22!j± m^lsj 

In our model computation, A; jS will become [28] 

A, s (T) = y^> (25) 

(O (o) - (t) (°») (°) 

In the right panel of Fig. 2, we present the comparison of A; jS obtained in PQMVT with LQCD [43,45,46] for 
both the parameter sets. We also compare with results from PQM. A similar analysis made in [28] concluded that 
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Figure 2: Left: Polyakov condensate (<&) obtained for ModelWB (solid blue) and ModelHotQCD (solid red) with the 
Polyakov potential Up \ y .yM- The dashed blue curve is ($) obtained in ModelWB with Up i y as the Polyakov potential. 
For comparison we have also plotted (<£>) obtained in PQM (solid black) with ModelHotQCD parameter set and Up \ y 
Polyakov potential. LQCD data from WB collaboration in blue (continuum estimate) and green (N T — 8) [45] and 
from HotQCD group in red (N T — 8 with asqtad action) and orange (N T = 8 with p4 action) [43] are also shown. 
Right: Plots of A; jS in ModelWB (blue) and ModelHotQCD (red) as obtained in PQMVT are shown. The black curve 
is the model prediction in PQM with ModelHotQCD parameter set. All the model predictions are obtained with the 
Polyakov potential Wp i y _vM ■ For comparison LQCD data are also shown: WB continuum estimate in blue [45] and 
N T = 10 [46] data in green are shown, the HotQCD data for N T = 8 with p4 (orange) and asqtad (red) actions [43] 
are also shown. 



the transition of A;. s is sharper in PQM than that obtained in lattice. Here we find that the inclusion of the vacuum 
term in PQMVT makes the drop in A/ iS across the crossover region much more gentle. This leads to better agreement 
with lattice data. Thus the crossover transition becomes smoother with the addition of the vacuum term. We will see 
that the smooth variation of the condensates across the transition temperature also result in smooth variation of all 
thermodynamic quantities and has important consequence on the phase diagram as well. 



3.1 p, e, A, s, c v , c 2 s 

While on the lattice, the usual convention is to compute the trace of the energy momentum tensor [43] and extract 
the rest of the thermodynamic quantities from it, in model calculations a different approach is followed. Here it is 
natural to first compute the grand potential f2 with the values of the condensates obtained by solving (20). From 
all the thermodynamic quantities can be easily obtained. The pressure of the system is given by 

p(T,fJ, x ,fly) = -CI (T, fJ, x , fly) (26) 

with the vacuum normalization p (0, 0, 0) = 0. The energy density e is given by 

*0(fi/T)_ rrf^l 



e = -T<^ = -T{-) + n (27) 

while the entropy density s is now easily obtained from the thermodynamic relation at zero chemical potential 

e = -p + Ts (28) 

Using (26) and (27), we can find the conformal symmetry breaking measure A given by the trace of the energy 
momentum tensor MM (T). 

= ^ (29) 

In Fig. (3), we have plotted these thermodynamic quantities as obtained in PQMVT for both the parameter sets, 
ModelWB and ModelHotQCD with Up \ y -vu as the Polyakov loop potential. We have also plotted the PQM results 
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Figure 3: Plots of p, e, s and A in PQMVT with ModclWB (blue solid) and ModclHotQCD (red solid) parameter sets. 
The corresponding black curves are the PQM predictions with ModclWB parameter set. These are obtained with the 
Polyakov potential Up \ y .vM- For comparison, model predictions with the Polyakov potential Up \ y and ModelWB 
parameter set are also shown in dashed blue in each case. Here T c = Tx ^ Td . The LQCD data from WB are plotted 
in green and blue while the HotQCD data arc plotted in orange and red. In case of p, the continuum estimate in 
blue and N T = 10 WB data in green [44] are shown. Also shown are the N T = 8 (red) [43] and N T — 6 (orange) [42] 
HotQCD data. For both e and s, LQCD data from WB are shown for N T = 10 (blue) and N T = 8 (green) [44] and 
HotQCD N T = 8 [43] (red) and N T = 6 [42] data are also shown. For A, WB data of continuum estimate [44] in blue 
and N T = 10 [46] in green are shown while HotQCD N T = 8 data with p4 (orange) and asqtad (red) [43] actions are 
also shown. In all the curves the high T SB limit of an ideal gas of massless fermions and gluons as in (30) is indicated 
by a thick black line while the dashed black line indicates the SB limit of an ideal gas of fermions only as in (30) 
without the first term. 
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with Wp iy_vM and PQMVT results with Up \ y for comparison. We find good qualitative agreement overall with LQCD. 
The smooth rise of p, e and s with T for model as well as LQCD indicate the crossover nature of the phase transition 
that takes place at [i = [28,38-41]. At low T, we see that $ and $ ~ and so it is clear from (17) that most of 
the contribution to thermodynamic quantities like p, e and s come from the three-quark "hadron" states. This is how 
statistical confinement works generically in these model studies [48]. As is evident from Fig. (1), high values of a x and 
a y result in large mass for these effective degrees of freedom yielding suppressed values for all these thermodynamic 
quantities at lower temperatures. At larger temperatures, the melting of the quark condensates across T x (see Fig. (1)) 
result in lighter effective degrees of freedom. Also, $ and $ — > 1 resulting in the release of color states as the single 
quark states become the dominant degrees of freedom. This is manifested by all the three thermodynamic quantities 
that steadily rise across T x and saturate to a much higher value close to that of an ideal gas of Nf massless fermions 
and (N 2 — l) gluons whose pressure psb is given by 

P ^ = W-l)i+WfW - (30) 

There is good agreement between ModelWB predictions for p, e and s with the corresponding WB N T — 8, 10 and 
continuum estimate [44] data in the temperature range w 0.8 — 2 T x . At lower temperatures the agreement between 
model predictions and LQCD is not good which could be due to the mean field treatment of the mesons. Comparison 
of LQCD data with model predictions after taking into account mesonic fluctuations [49] should be carried out in 
order to throw more light on this issue. At higher temperatures the model description may not match well with 
LQCD as the contribution from the transverse gluons become significant and their physics may not be captured by 
the Polyakov loop $. The ModclHotQCD prediction although qualitatively in agreement with LQCD, lies slightly 
below the HotQCD data for all temperatures. However, it should be noted that the N T = 8 data [43] is consistently 
lower than that of N T = 6 [42] and so we expect that the agreement will improve when the continuum estimates will 
be available. 

The model predictions for the conformal symmetry breaking measure A has also been plotted and compared with 
LQCD. In this case we have plotted A against T c = Tx ^ Td for better agreement between model and LQCD. While 
for lower and higher temperatures, HotQCD [43] and WB [44,46] data look reasonably consistent, around 1 — 1.5 T c 
although both exhibit a peak, the height of the peak in case of HotQCD is almost 50% greater than that of WB. 
Moreover, the N T = 8 HotQCD data with asqtad action is significantly lower than that of p4 action at the peak. Thus 
the LQCD prediction in this temperature regime is highly sensitive to the lattice volume as well as the lattice action 
implemented. As mentioned earlier, we have chosen the value of k in (8) such that the model prediction agrees well 
with LQCD. In general k could be taken as a function of T but in this work we have taken a constant value for k. 

The bumpy structures in (<r y ) and ($) that we found around T x in Fig. 1 in the case of PQM shows up in the 
thermodynamic quantities also. Although the pressure looks smooth across the transition region, others like the e, s 
and A which are functions of |p exhibit bumps around T x . As observed in Fig. (1), addition of the vacuum term in 
PQMVT results in a smooth behaviour of the condensates in the transition region which in turn smoothens all the 
thermodynamic quantities in Fig. (3). 

Another important quantity from the point of view of hydrodynamical investigations of relativistic heavy ion 
collisions is the isentropic speed of sound c s given by 



dp 


dp 




/ de 


s 


de 


s = df 


vl 


df 


v c v 



where cy is the specific heat capacity at constant volume obtained as follows 

d 2 n 



de 

CV= df 



-T 

g T 2 



v 

c 2 is also closely related to the equation of state parameter p/e [43] 



(31) 



(32) 



c 2 = ^ = el (£) + (£) (33) 



de de 



In Fig. 4 we have plotted the model prediction for c 2 and p/e against T and also e. In the latter case the uncertainties 
due to T x go away. In Fig. 5 we have plotted the model prediction for c v . The bumpy structures seen in Fig. 3 for 
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Figure 4: Plots of c 2 s (solid) and p/e (dashed) are shown in PQMVT with ModelWB (blue) and ModelHotQCD (red) 
parameter sets. The corresponding black curves are the PQM predictions with ModelWB parameter set. These are 
obtained with the Polyakov potential £/p i y _vM- Also shown are the plots with Up i y with ModelWB parameter set in 
brown for comparison. WB data of c 2 s for N T = 10 in green and the continuum estimate in blue [44] are shown. Also 
shown are the HotQCD N T = 8 data [43] of p/e in red. In all the curves the SB limit is indicated by a thick black line. 




Figure 5: Plots of c v in PQMVT with ModelWB (blue) and ModelHotQCD (red) are shown. The black curve is the 
PQM prediction with ModelWB parameter set. These are obtained with the Polyakov potential Wp i y _vM- Also shown 
is the plot with Up \ y with ModelWB parameter set (blue dashed) for comparison. The high T SB limit of an ideal 
gas of of massless fermions and gluons as in (30) is indicated by a thick black line while the dashed black line indicates 
the SB limit of an ideal gas of fermions only as in (30) without the first term. 
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PQM becomes even more prominent in c 2 s and cy which depend on cy exhibits a sharp spike like structure at 

T x which results in a sharp dip in c 2 s at the same temperature. The PQMVT plots are lot smoother and compare well 
with LQCD qualitatively. Around T x , WB continuum estimate for <? s [44] show a minimum with a value around 0.12 
while ModelWB also exhibits a minimum in the same region with a minimum value ~ 0.08. In the high T limit, both 
c 2 and p/e for model and LQCD are found to approach the SB limit 1/3. In case of cy the high temperature limit is 
dependent on the Polyakov potential chosen. As seen in Fig. (5), the model prediction with Up i y saturates to the SB 
limit by 2T X whereas the model prediction with Wp i y _vM saturates to a lower value ~ 73% of the SB limit. Similarly 
even in the case of p, e and s in the high temperature we find the model predictions with £/p i y _vM always lower than 
those with Up \ y . 

We end this subsection with an interesting observation. We find that the model prediction for the thermodynamic 
quantities like p, e, s and cy with ModelWB parameter set and Wp i y -VM Polyakov potential for high temperatures 
saturate close to the SB limit of an ideal gas with only massless fermions and no gluons which is the psb from (30) 
after removing the first term due to gluons. In [28] similar results were reported when the Fukushima ansatz [50] 
for the Polyakov potential was used. Thus we find that the high temperature behaviour for model predictions with 
ModelWB parameter set and Wp i y -VM is similar to that of the Fukushima potential [50] which is based on a strong 
coupling analysis. 



3.2 Quark Number Susceptibilities 

Quark number susceptibilities (QNS) are defined as derivatives of pressure with respect to quark chemical potentials 

Xuk = ~ ^ ~ b (34) 

Linear combinations of these can be related to fluctuations and correlations of conserved charges like baryon number, 
electric charge and strangeness which have been proposed as possible observables to signal the QCD phase transi- 
tion [54,55]. These have been measured in LQCD [43,45,47,51-53]. Here we show results on QNS upto second order. 
A detailed analysis of the fluctuations and correlations of conserved charges in PQMVT and their comparison with 
LQCD will be done elsewhere [56]. 

Fig. (6) shows the diagonal susceptibilities x% an d X21 their ratio X2/X2 an d the off-diagonal xT '■ All susceptibilities 
vanish at low temperature since the condensates have large values which implies large masses of the relevant degrees 
of freedom that restricts the capacity of fluctuations. At high temperature where the degres of freedom have small 
masses and temperature is the only relevant scale, the system begins to behave like an ideal gas of massless quarks. 
The susceptibilities thus approach the corresponding SB limit of ideal gas of three quarks. 

The diagonal susceptibilities both show a monotonic increase with temperature, with the maximium rate of increase 
observed around the crossover region and then saturating at higher temperatures to the SB limit. Both X2 an(1 X2 
saturate to about 80 — 90% of the SB limit at T <~ 2.4T X . We also observe that the light quark susceptibility rises 
much faster than that of the strange quark both in model as well as LQCD. xl begins to saturate well beyond the 
cross-over region as is expected from the slow decrease of the value of the strange condensate. For both X2 an d X2 
the N T = 8 data with asqtad action for HotQCD [43] lie consistently above ModclHotQCD. On the other hand, there 
is much better agreement between WB data and ModelWB predictions. It is interesting to observe that the recent 
N T = 8 susceptibility data with hisq action from HotQCD [47] which has also been plotted in Fig. (6) agree much 
better with the WB data. 

It is useful to consider the ratio of the two susceptibilities, X2/X2 > m order to avoid unknown normalization factors 
that maybe present in lattice data [16,52]. In the unbroken phase there is excellent agreement between the models 
and lattice data, but they do not match as well below the crossover temperature. This maybe attributed to absence of 
mesonic fluctuations in a mean field approximation calculation [28]. However, this ratio in the low temperature region 
is also highly dependent on the form of the Polyakov potential (see Fig. (6 b) of [28]). In our case the potential we 
use improves marginally the agreement of the model to lattice predictions by giving smaller values to this ratio in the 
low temperature region. It should also be noted that the lattice data in this region have large errors and comparisons 
of the model to lattice are uncertain. 

We also show in Fig. (6) the off-diagonal susceptibility xT S which represents the correlation between the light and 
heavy quarks. It vanishes in the SB limit, as is expected for a non-interacting ideal gas, and peaks in the crossover 
region. Qualitatively the X2 S predictions of the model agree with lattice, showing the same essential features, however, 
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Figure 6: QNS as obtained in PQMVT are plotted for ModelWB (blue) and ModclHotQCD (red) with the Polyakov 
potential W Po i y _ V M- For both X % and xl, WB N T = 12 data [51] (green) and HotQCD N T — 8 [43] (red), N T = 6 [52] 
(orange) and N T = 8 with hisq action [47] (brown) are shown for comparison. Also in case of xIj WB continuum 
estimate [45] is shown in blue. The HotQCD data for X2/X2 m re d [43] and the WB continuum estimate for xY [53] 
in blue are also compared with the model predictions. The high T SB limit is indicated by a thick black line. 
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Figure 7: Phase diagrams for PQM (black) in ModelWB and PQMVT in ModelWB (blue) and ModelHotQCD (red) 
with the Polyakov potential Wp i y -VM- The dashed lines indicate crossover, the solid lines indicate first order phase 
transition and the filled circles indicate CEP. The phase diagram for ra„ — 600 MeV in the case of PQMVT with the 
Polyakov potential hip \ y is also shown (green) . 

lattice data is consistently larger over the entire temperature range. This maybe due to unknown normalization factors 
of lattice and perhaps a better quantity to compare would be normalized correlations [56] . 

3.3 Phase Diagram 

The divergence of %2 can be used to locate the CEP [57] where a first order phase transition line ends in a second 
order phase transition point. We obtain the phase diagram in the case of symmetric quark matter (fi x = \x y = fi). In 
Fig. (7) we have presented the phase diagram for PQM and PQMVT. The chiral phase boundary has been plotted 
and the CEP located in all the cases. The smoothening of the thermodynamic quantities in the transition region as 
seen for /x = persists even at non-zero chemical potentials. This results in pushing the CEP to a higher value of ii in 
case of PQMVT as compared to PQM. With m a = 400 MeV, in PQM with ModelWB using the W Poly . V M potential 
the CEP is found at (jj,,T) = (80, 177.5) MeV which gets shifted to (283,81) MeV in case of PQMVT of ModelWB. 
We have also obtained the phase diagram in PQMVT with m a = 600 MeV and U-p \ y potential. In this case the phase 
diagram has no CEP and there is only crossover transition on the entire (/x — T) plane as in the case of 2 flavors [26] . 

We have also plotted the phase boundary of PQMVT in ModelHotQCD (red) with the Wp i y _vM potential. In this 
case T x is lower than that of ModelWB at fj, = that result in a lower phase boundary in the (T, /x) plane as compared 
to the case with ModelWB parameter set. This is the effect of choosing a smaller value of Tq. For large values of fx 
there is not much difference between the two phase boundaries. This is essentially because in this region the Polyakov 
potential plays a less important role as compared to Vt q g. The phase boundary for ModelWB with Up \ y does not differ 
much from that of ModelWB with the Wp i y _vM potential implying that in this case the VanderMonde term does not 
play much of a role in the determination of the chiral phase boundary. 

4 Conclusion 

In this paper we have studied the effect of including the commonly neglected fermionic vacuum fluctuations to the 
(2 + 1) PQM model. The conventional PQM model suffers from a rapid phase transition contrary to what is found 
through lattice simulations [28] . This could be due to the use of a Polyakov potential which carries with it a remnant 
first order phase transition of the pure glue theory. Another possible reason could be the sudden release of quark 
degrees of freedom in PQM with increase in temperature. Addition of the vacuum term in PQMVT addresses the latter. 



13 



This tames the rapid transitions that we see in the PQM model and significantly improves the model's agreement to 
lattice data. 

We work with two Polyakov loop potentials: Up \ y and Wp i y -VM given in (7) and (8) respectively. The model 
parameters are fixed by inputs from experiments as well as LQCD data. In this work, we have compared model 
predictions to LQCD data of HotQCD and WB groups. Depending on the lattice group whose data we choose to fix 
the model parameters, we arrive at the two sets of model parameters namely ModelWB and ModclHotQCD given 
in Table 1. We then investigated the bulk thermodynamic properties of the system and the influence of the vacuum 
term in PQMVT. In PQM, A;. s which is a suitably defined order parameter for chiral phase transition is found to 
undergo a much rapid fall across the transition temperatures compared to LQCD [28]. The inclusion of the fermionic 
vacuum fluctuations in the model makes the transition much more gentle. The vacuum term essentially washes away 
all the jagged and bumpy structures that are there in the quark codensates, Polyakov loop and all thermodynamic 
quantities in PQM. We also computed the QNS upto second order. In all cases very good qualitative agreement 
between LQCD and model is found. The model computations suggest that the high temperature behaviour of the 
thermodynamic quantities with Wp i y -VM Polyakov potential and ModelWB parameter set is similar to that of the 
Fukushima potential [50] for the Polyakov loop. We further investigated the role of the vacuum term on the phase 
diagram. The smoothening effect of the vacuum term persists even at non zero fi. With m a = 600 MeV, we find that 
the vacuum term washes away the critical behaviour completely from the (fi, T) plane. In case of m a = 400 MeV, 
addition of the vacuum term pushes the CEP to larger chemical potential. 
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6 Appendix 

A Model Parameters 

The model parameters are determined by establishing the vacuum properties. In analogy to [15,32] the values of the 
condensates are determined from the pion and kaon decay constants by means of the partially conserved axial-vector 
current relation (PC AC). In the strange-non-strange basis they are given by 

o x = U I &y = ~y= ( 2 /if - U) ■ ( 35 ) 

The remaining model parameters (to 2 , Ai, A2, c, h x and h y ) are fixed by the following inputs from experiments: pion 
mass TOtt and pion decay constant f v , kaon mass vhk and kaon decay constant Jk, average squared mass of 77 and 
rf mesons, (m 2 + to 2 ,) and sigma mass m a . In this appendix we shall derive the equations to determine the above 
parameters in PQMVT. 

Using (10) and (14), we rewrite fi (T 7 fi X7 fi y ) as 

qPQMVT (T) ^ ^ = PQM (T) ^ ^ + ^ (36) 

where 

fi p Q M (T, fi Xlt i y ) = U(a x ,a y ) + U P ($, $, T) + fig (37) 
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Mesons 



Expression for masses 



2 + Ar(.x 2 + y 2 ) + 4f (x 2 - V2xy + 2y 2 ) - fx - gg£ (fpgf) (x 3 X + 2^Y) 
• 2 + X 1 (x 2 + y 2 ) + % (x 2 + y 2 ) + §(2x + ^y) - ^ (x 2 X + y 2 Y) 
2 + AiCx 2 + y 2 ) + f (x 2 + 4y 2 ) - §(4x - V2y) - ^ (x 2 X + Ay 2 Y) 

'p. 00 + m p.88 

,2 + ^ (7 . T 2 + A ^2 xy + 5?y2) + A2(x 2 + y2) _ v|c ( 72x + y) - ^ (3 (x 2 X + y 2 Y) + 4 (x 2 + y 2 )) 
,2 + Ax(5. T 2 _ ^y + 7y 2 ) + A 2 (4 + 2y 2 ) + ^(^x - § ) - ^ (| (x 2 X + 4y 2 y) + 2 (x 2 + 4y 2 )) 

9 A 



to k 

2 

m p,00 
2 

TO p,88 

(to 2 4 

9 

TO 



s.00 



TO 

m 

TO 
TO 



^(^x 2 - xy - vV) + V2A 2 (^ - y 2 ) + ^(x - ^2y) -$£r Q (x 2 X - 2y 2 F) + § (x 2 - 2y 2 )) 
TO loo cos2 0s + m-s,88 sin2 ^ + 2m 2 08 sin6> s cos6> s , where tan (26 8 ) = m i 2 ™^™i ^ 



Table 2: The expressions of squared meson masses in the vacuum required to determine the model parameters, x 
denotes u Xl y denotes a y , X denotes (l + 4 log (^f )) and Y denotes ^1 + 4 log (jj^j^- 



Now the meson masses are determined by the curvature of f2 at the global minimum 

d 2 n 



m 



a,ab 



5 2^PQM 



+ 



d 2 m 



QQ 



PQM 2 . v 2 
m *ab + m a ,ab 



d^a,ad^ a ,b 



are available in the literature [32]. Here we compute rrC a a 2 . From (19) and (38) we get 



L a,ab 



~8tt 2 ^ 



f L 



2 log 



dm 2 



dm 2 



^£ +2 TO 2 l0g(^ 



\ \ m f 

> ) d£ a , a d£ a ,b 



(38) 



(39) 



The expressions for (^ g™ f ^ ^ g"/ ) anc ^ eg 8 ™| b are already available in [32]. In Table (2), we present the expressions 
of to 2 ab for the different mesons as obtained in PQMVT at T = = 0. 
From the expressions of to 2 , m 2 K and (to 2 + to 2 ,) we find A2 and c 



where 



A 2 = 



c = 



m 



3 (2f K - U) m' 2 - (2f K + U) m' 2 - 2 (to 2 + to 2 ,)' (f K - f w ) 



' 1 '2 

!k — f-rr 



(m+8fK(f K -u))(f K -u) 

A 2 (2f K - U) 



(40) 



m 



"1 



TO, 



64tt 2 

2 1 _2 



(41) 



96tt 2 v y ; 192tt 2 

Here x denotes er x , y denotes cr y , X denotes (l + 4 log (^f )) and Y denotes ^1 + 4 log {^j^^j- Having obtained A 
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and c, we use the expression of m 2 . to express m 2 in terms of Ai 

- 2 = m 2 - ^fl + | (2/k - /tt) - Ai [fl + (2/ " 2 ' A)2 ) (42) 

Finally, Ai may be obtained from the expression of m 2 while the explicit symmtery breaking parameters h x and /i^ 
are obtained from the stationarity conditions (20) 



V2A 



(44) 
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Thus the final expression of the parameters are 
h x 

hy 

A 2 (A) 



m 



m t,oo 



''s,08 



tan (26 a ) = 







2fKm 2 K 


-fa 


V2 


(^.9 4 Log 


~gU~ 




2A 



ff 4 Log 



.f*(-8f 2 K + 8.f K .U-3f*)N c + 
g(2f K -UY 



2A 



(-2f K + U) 2 (8/| - + 3/ x 2 ) iV c + 



Sk (8g 4 f K N c - 16g 4 f K UN c - f v (32tt 2 (3m 2 , + m 2 - 2 (m 2 + m 2 ,)) + 3g 4 f 2 N c ) + 
Ik (-64^ 2 (-3m 2 , + m 2 + m 2 + m 2 ,) + Ug 4 f 2 N c ))) / (32v 2 f K (8/| - 16/|/ OT + ll/„/ 2 - 3/ 3 )) ; 
2 (2/g (m 2 , + m 2 - m 2 - m 2 ,) + / ff (-2m 2 + m 2 + m 2 rj ,)) _ 



8/ 2 



A' 



8/a-/ w + 3/ 2 



-/Log 



f 2 (-2f K + U) 2 (8f 2 K 



K U + 3/ 2 ) N c 

2 



32tt 2 /„ (32/^Ai - 96/i/^Ai + / 3 (-3m 2 , + m 2 + 4 (m 2 + m 



9/^Ai) 



voo 



COS 



= 3 ff 4 Log 



(2/k ^ U) 



tficU (3m 2 , - 3m 2 - 4 (m 2 + m 2 ,) + 26/ 2 Ai) + Af% (-2 (-m 2 , + m 2 + m 2 + m 2 ,) + 33/ 2 Ai) + 
/*/ 2 (10m 2 , - 8m 2 - 12 (m 2 + m 2 ,) + 45/ 2 Ax))) / (64^ 2 / i , (8/| - 16/1/, + Uf K f 2 - 3/ 3 )) ; 
2 S + 2m 2 08 sin S cos # s , where 

f 2 (-2f K + U) 2 (8/| - 8f K .U + 3/ 2 ) N c - 

2f K (32/| ( 5 4 iV c - 16^%) + f K U (-96.g 4 7V c + 51277%) + 

f K f 2 (16tt 2 (-26m 2 , - 9m 2 + m 2 + m 2 ,) + f 2 (45g 4 N c + 16w 2 X 1 )) - 

8f 2 K .U (IQtt 2 (-5m 2 , - 3m 2 + m 2 + m 2 ,) + / 2 (l3 5 4 7V c + 32^)) + 



fA-^n 2 



-9m 



K 



2m 2 + m 2 + m 2 ,) + f 2 (-9g 4 N c + 48tt 2 Ai)) + 



4/| (32tt 2 (-4m 2 , - m 2 + m 2 + m 2 ,) + / 2 (33.g 4 iV c + 4877%)))) / (96^ 2 /a- (8/| - 16/|/,+ 
ll/*/ 2 - 3/ 3 )) ; 



= 3g 4 Log 



fl (-2f K + U) 2 (8/| - 8f K .U + 3/ 2 ) Nc- 



Va (-64/4 / w _ 64 ^2 Ai ^ + 64/ 5 ( g i Nc _ Wn 2 X ^ + 

Ml (16^ 2 (-58m 2 , - 39m 2 + 32 (m 2 + m 2 ,)) + / 2 (57 5 4 Ar c - 217671%)) - 

32 f 2 K U (4tt 2 (-13m 2 , - 3m 2 + 5 (m 2 + m 2 ,)) + / 2 (5 ff 4 iV c - 164tt 2 A 1 )) + 

16/| (8tt 2 (-8m 2 , + m 2 + 2 (m 2 + m 2 ,)) + 3/ 2 (5.g 4 Ar c - 13677%)) + 

fl (16tt 2 (9m 2 , + 16m 2 - 8 (m 2 + m 2 ,)) + 

fl (-9g 4 N c + 384tt%)))) / (96^/k (8/| - 16/ 2 -/, + llf K f 2 - 3/ 3 )) ; 
(8/i (g 4 N c - 16t7%) - 16/|/ T (g 4 N c - 12tt 2 Ai) + 8tt 2 / 2 (2m 2 - m 2 - m 2 , + 3/ 2 Ai) + 
f K (32tt 2 (-4m 2 , + 2m 2 + m 2 + m 2 ,) + / 2 (ll ff 4 iV c - 487r 2 A!)) - 

fxU (8tt 2 (-7m 2 , + m 2 + 3 (m 2 + m 2 ,)) + / 2 (3.g 4 iV c + 40^))) / (e^Tr 2 (8/| - 8/^/. + 3/; 



2TO ^08 



(45) 



s,00 



''s.88 



Hence only A2 depends on A. It is straightforward but tedious to check that this dependence cancels neatly with that 
of in (10) to yield a A independent £1 
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